Positive stochastic volatility simulation

We present a positivity preserving numerical scheme for the pathwise solution of nonlinear stochastic differential equations driven by a multi-dimensional Wiener process and governed by non-commutative linear and non-Lipschitz vector fields. This strong order one scheme uses: (i) Strang exponential splitting, an approximation that decomposes the stochastic flow separately into the drift flow, and the pure diffusion flow governed by the diffusion vector fields; (ii) an implicit Euler method to approximate the drift flow; and (iii) an implicit Milstein method to approximate the pure diffusion flow. The separate approximations for the drift and pure diffusion flows preserve positivity. Therefore the Strang exponential splitting approximation does also. We demonstrate the efficacy of our method by applying it to the Heston model and a variance curve model, and compare it against well-established positivity preserving schemes.